How frequently revisited sites differ between individuals and the sex of the spotted nutcracker in Davos

Project for Patterns and Trends in Environmental Data - 2025

pacman::p_load(
  adehabitatHR,
  cluster,
  dplyr,
  forcats,
  ggplot2,
  ggpubr,
  ggridges,
  knitr,
  lubridate,
  multcompView,
  patchwork,
  pacman,
  pracma,
  purrr,
  readr,
  recurse,
  sf,
  terra,
  tidyr,
  viridis,
  wordcountadding
)


set.seed(123)

theme_set(
  theme_bw()
)

scale_lv95_axes <- function() {
  list(
    scale_x_continuous(
      limits = c(2781000, 2789000),
      minor_breaks = c(2781000, 2783000, 2785000, 2787000, 2789000),
      breaks = c(2782000, 2784000, 2786000, 2788000)
    ),
    scale_y_continuous(
      limits = c(1185000, 1195000),
      minor_breaks = c(1185000, 1187000, 1189000, 1191000, 1193000, 1195000),
      breaks = c(1186000, 1188000, 1190000, 1192000, 1194000)
    ),
    coord_sf(
      datum = st_crs(2056),
      expand = F
    ),
    theme(
      axis.text.x = element_text(angle = 45, hjust = 1),
    )
  )
}

scale_animal_id <- function() {
  list(
    scale_color_brewer(
      palette = "Spectral",
      name = "animal id",
    ),
    scale_fill_brewer(
      palette = "Spectral",
      name = "animal id",
    )
  )
}

scale_animal_sex <- function() {
  list(
    scale_fill_brewer(
      palette = "Spectral",
      name = "animal sex",
      labels = c("female","male")
    ),
    scale_color_brewer(
      palette = "Spectral",
      name = "animal sex",
      labels = c("female","male")
    )
  )
}

background <- png::readPNG("data/swissimage2.5m_latest.png")   

spotted_nutcrackers_ref_data <- read_delim("data/reference_data_spotted_nutcrackers.csv") |>
  mutate(across(where(is.character), as.factor)) |> 
  mutate(
    tag_id = factor(tag_id)
  )

spotted_nutcrackers_gps <- read_delim("data/gps_spotted_nutcrackers.csv") |>
  mutate(
    tag_local_identifier = factor(tag_local_identifier),
    timestamp = as.POSIXct(timestamp, format = "%d/%m/%Y %H:%M"),
    across(where(is.character), as.factor))

spotted_nutcrackers <- left_join(
  spotted_nutcrackers_gps,
  spotted_nutcrackers_ref_data,
  join_by(individual_local_identifier == animal_id)
) |> 
  st_as_sf(coords = c("location_long", "location_lat"),
           crs = 4326) |>  # Setzt ursprüngliches CRS
  st_transform(crs = 2056) |>  # Transformiert ins LV95
  mutate(CH_X = st_coordinates(geometry)[, 1], # Extrahiert neue X-Koordinaten
         CH_Y = st_coordinates(geometry)[, 2]) |>
  st_drop_geometry()

Background and Research Goals (simon)

Nutcrackers (Nucifraga caryocatactes) play a crucial role in forest regeneration by dispersing seeds. Analyzing their movement data allows us to identify frequently revisited sites, which can provide insights into harvesting and caching behavior as well as dispersal distances. Understanding these movement patterns and their ecological implications is essential for forest conservation and management (Sorensen et al. 2022; Graf et al. 2024).

1. How do frequently revisited sites vary among individuals, and can distinct clusters of these sites be identified?
2. Do revisited sites of individuals differ between males and females in terms of distance from each other?

Material and Methods

Study Area

Both studies (Sorensen et al. 2022; Graf et al. 2024) tracked spotted nutcracker movements in 2017 and 2018 within a ~15 km² area around Davos in the eastern Swiss Alps (Figure 1). In this region, Swiss stone pine (Pinus cembra) forms the upper treeline between 1850 and 2200 m a.s.l., with the highest densities at mid-elevations. Valley bottoms are dominated by Norway spruce (Picea abies) and European larch (Larix decidua), where Swiss stone pine is scarce.

Figure 1: Study Area around Davos

Data

The data used consists of two datasets: one containing GPS localizations with x/y coordinates, timestamps, and individual identifiers, and a reference dataset with additional information on each individual, such as sex, weight, and wingspan. In total, data from 31 individuals were recorded, tracked over varying time periods and at different temporal resolutions, with localization intervals ranging from 15 to 60 minutes. We only used data collected in the months of august and september in 2017 and 2018. Individuals were tracked between 10:00 and 20:00 with a 15 minute frequency. For the identification of the freqeuntly revisited sites, we used all gps localizations, for the analysis of the stepwidth only those with a time lag of 15 minutes to the previous localization.

spotted_nutcrackers_filter <- spotted_nutcrackers |> 
  group_by(animal_ring_id) |> 
  summarise(
    count = n(),
    start =first(timestamp),
    last=last(timestamp),
    duration = last - start 
  ) |> 
  filter(
    duration > 5,
    start |> year() %in% c(2017,2018),
    last |> year() %in% c(2017,2018),
    !animal_ring_id %in% c("K96915","K96934")
  ) |> 
  pull(animal_ring_id)

spotted_nutcrackers_filtered <- spotted_nutcrackers |> 
  filter(animal_ring_id %in% spotted_nutcrackers_filter) 

spotted_nutcrackers_filtered |> 
  dplyr::select(animal_ring_id, timestamp, CH_X, CH_Y, animal_sex) |> 
  head() |> 
  kable()

Table 1: excerpt of the data used

animal_ring_id	timestamp	CH_X	CH_Y	animal_sex
K125864	2018-08-10 16:15:00	2788349	1186097	f
K125864	2018-08-10 16:45:00	2788316	1186106	f
K125864	2018-08-10 17:15:00	2788242	1186158	f
K125864	2018-08-10 17:45:00	2787367	1186915	f
K125864	2018-08-11 10:00:00	2787844	1186519	f
K125864	2018-08-11 11:00:00	2787889	1186550	f

Method

To identify highly revisited sites, we applied the workflow described by Sorensen et al. (2022), conducting a recursion analysis for each bird’s GPS locations using the R package recurse (Bracis, Bildstein, and Mueller 2018). Each GPS point was treated as a potential revisitation site by defining a 250 m radius around it. We then counted the number of revisits, i.e., instances where a bird returned to the same site. The density distribution of revisits often exhibited a multi- or bimodal pattern. The lowest local minimum was used as a threshold to distinguish points with few or no revisits from frequently revisited points.

find_local_minimum <- function(values) {
  d <- density(values)
  dy <- d$y
  dx <- d$x
  minima_idx <- which(diff(sign(diff(dy))) == 2) + 1
  if (length(minima_idx) == 0) return(NA)
  minima_x <- dx[minima_idx]

  lowest_minimum <- minima_x |> min()
  
  return(lowest_minimum)
}

spotted_nutcrackers_recurse <- spotted_nutcrackers_filtered |>
  group_by(animal_ring_id) |>
  group_map(~ {
    rec <- .x |> dplyr::select(CH_X, CH_Y, timestamp, individual_local_identifier) |> 
      as.data.frame() |> 
      getRecursions(radius = 250, timeunits = "mins")
    rec$revisitStats
  }) |> 
  map_dfr(~ .x) |> 
  group_by(id, coordIdx, x, y) |> 
  summarise(revisits = n(), .groups = "drop") |> 
  group_by(id) |> 
  mutate(
    threshold = find_local_minimum(revisits),
    mean = mean(revisits),
    median = median(revisits),
    max_revisits = max(revisits),
    freq_revisited = revisits >= threshold
  ) |> 
  ungroup() |> 
  left_join(
    spotted_nutcrackers_ref_data,
    by = join_by(id == animal_id)
  )

The frequently revisited points then underwent a spatial fuzzy cluster analysis. For each individual they were grouped into two distinct clusters using the R package cluster (Maechler et al. 2024). If the two centorides of the two clusters were closer than 250 m together, they were combined into a single cluster. Extracting the 95% minimum convex polygon for each cluster using the adehabitatHR package (Calenge 2024) resulted in one or two frequently revisited sites for each bird.

spotted_nutcrackers_points <- spotted_nutcrackers_recurse |> 
  st_as_sf(coords = c("x", "y"), crs = 2056)

spotted_nutcrackers_clustered_points <- spotted_nutcrackers_recurse |> 
  filter( freq_revisited) |> 
  dplyr::select(id, x, y, animal_sex,  freq_revisited) |> 
  group_by(id) |> 
  group_modify(~ {
    coords <- dplyr::select(.x, x, y)

    if (nrow(coords) < 2) {
      .x$cluster <- NA
    } else {
      coords_sf <- st_as_sf(coords, coords = c("x", "y"), crs = 2056)
      
      pam_result <- pam(st_coordinates(coords_sf), k = 2)
      .x$cluster <- pam_result$clustering

      cluster_centroids <- coords_sf |>
        mutate(cluster = .x$cluster) |>
        group_by(cluster) |>
        summarise(geometry = st_centroid(st_union(geometry))) |>
        st_as_sf()

      dist_between_clusters <- st_distance(cluster_centroids[1, ], cluster_centroids[2, ])

      if (as.numeric(dist_between_clusters) <= 250) {
        .x$cluster <- 1
      }
    }

    return(.x)
  }) |> 
  ungroup() |> 
  mutate(id_cluster = paste(id, cluster, sep = "_")) |> 
  group_by(id_cluster) |> 
  filter(n() >= 5) |> 
  ungroup() |> 
  st_as_sf(coords = c("x", "y"), crs = 2056) |> 
  as("Spatial")

spotted_nutcrackers_clustered_mcp <- spotted_nutcrackers_clustered_points["id_cluster"] |> 
  mcp(percent = 95) |> 
  st_as_sf(crs = 2056) |> 
  separate(id, into = c("id", "cluster"), sep = "_") |> 
  arrange(area |> desc()) |> 
  left_join(
    spotted_nutcrackers_ref_data,
    join_by(id==animal_id)
  )

To quantify movement activity, we calculated the stepwidth (Euclidean distance in meters) between consecutive GPS locations for each individual. The stepwidth was computed using a custom function that calculates pairwise distances between subsequent locations. Additionally, the time lag between fixes was calculated in minutes, and only steps with a regular 15-minute interval were retained for analysis to ensure temporal consistency.

get_stepwidth = function(geometry, n = 1) {
  st_distance(lag(geometry, n = n), geometry, by_element = TRUE) |> as.numeric()
}

get_timelag = function(datetime, n = 1, units = "secs") {
  difftime(lead(datetime, n = n), datetime, units = units) |> as.numeric()
}

spotted_nutcrackers_stepwidth <- spotted_nutcrackers_filtered |> 
  st_as_sf(coords = c("CH_X", "CH_Y"), crs = 2056) |> 
  group_by(animal_ring_id) |>
  mutate(
    stepwidth = geometry |>  get_stepwidth(),
    timelag = timestamp |> get_timelag(units="mins"),
    speed_m_sec = stepwidth / (timelag * 60),
    speed_km_h = speed_m_sec * 3.6
  ) |> 
  ungroup() |> 
  dplyr::select(animal_ring_id, timestamp, stepwidth, timelag, animal_sex,animal_group_id) |> 
  filter(timelag==15) |>
  mutate(CH_X = st_coordinates(geometry)[, 1], 
         CH_Y = st_coordinates(geometry)[, 2]) |>
  st_drop_geometry()

To assess potential differences in movement behavior, we visualized stepwidth distributions with boxplots grouped by individuals and sex. We performed ANOVAs to test for significant differences in stepwidth between the individuals as well as males and females. Stepwidth therefore was transformed using the log10 function.

spotted_nutcrackers_aov <- aov(log10(stepwidth)~animal_ring_id, data=spotted_nutcrackers_stepwidth)

spotted_nutcrackers_tukey <- spotted_nutcrackers_aov |> TukeyHSD()
spotted_nutcrackers_aov_letters <- multcompLetters4(
  spotted_nutcrackers_aov,
  spotted_nutcrackers_tukey
)

spotted_nutcrackers_aov_label <- data.frame(
      animal_ring_id = names(spotted_nutcrackers_aov_letters$animal_ring_id$Letters),
      group = spotted_nutcrackers_aov_letters$animal_ring_id$Letters |> factor(levels = c("a","ab","abc","bc","c"), ordered=T))


spotted_nutcrackers_aov_sex <- aov(log10(stepwidth)~animal_sex, data=spotted_nutcrackers_stepwidth)

spotted_nutcrackers_points <- spotted_nutcrackers_points |> 
  left_join(
    spotted_nutcrackers_aov_label,
    join_by(animal_ring_id)
  ) |> 
  mutate(
    animal_ring_id = fct_reorder(animal_ring_id, group |> factor() |>  as.numeric())
  )

spotted_nutcrackers_stepwidth <- spotted_nutcrackers_stepwidth |> 
  left_join(
    spotted_nutcrackers_aov_label,
    join_by(animal_ring_id)
  ) |> 
  mutate(
    animal_ring_id = fct_reorder(animal_ring_id, group |> factor() |>  as.numeric())
  )

spotted_nutcrackers_clustered_mcp <- spotted_nutcrackers_clustered_mcp |> 
  left_join(
    spotted_nutcrackers_aov_label,
    join_by(animal_ring_id)
  ) |> 
  mutate(
    animal_ring_id = fct_reorder(animal_ring_id, group |> factor() |>  as.numeric())
  )

Results

The threshold value for distinguishing GPS localizations with infrequent revisits from those with frequent revisits ranges from 11.9 to 45.2, depending on the individual and their density curve of the number of revitis per GPS localization from the recursion analysis (Figure 2).

spotted_nutcrackers_points |>   
  ggplot(aes(revisits, fill=animal_sex)) +
  geom_density() + 
  facet_wrap(.~animal_ring_id, ncol=5,scales="free")+
  geom_vline(
    aes(xintercept = threshold),
    color = "blue", linewidth = 0.5
  )+
  labs(
    y="density"
  )+
  geom_text(
    aes(
      max_revisits/2,
      0,
      label = paste0("threshold: ",round(threshold, 1))
    ),
    color="blue",
    vjust = -0.5, 
    size = 3
  )+
  scale_animal_sex()

Figure 2: Density of revisitation of GPS licalizations per individual

The visualization of the frequently revisited GPS localization already shows a distinct pattern of frequently revisited sites (Figure 3). Clustering these points into two clusters and extracting the 95% minimum convex polygons for each individual results in two frequently revisited sites further than 250 meters apart for 9 out of 10 individuals. In one case, the two clusters are closer than 250 meters, so we identified only one frequently revisited location for this individual. All but one individual showed consistent revisits to one area in the southeast and another in the northeast or east.

p1 <- ggplot()+
  background_image(background)+
  geom_sf(
    aes(color=factor(animal_ring_id)),
    data=spotted_nutcrackers_points |> 
      filter(!freq_revisited),
    alpha=0.2
  )+
  geom_sf(
    aes(color=factor(animal_ring_id)), 
    data=spotted_nutcrackers_points |> 
      filter(freq_revisited)
  )+
  scale_animal_id() +
  scale_lv95_axes()+
  theme(legend.position = "none")

p2 <- ggplot() +
  background_image(background)+
  geom_sf(
    aes(color=factor(animal_ring_id)), 
    data=spotted_nutcrackers_points |> 
      filter(!freq_revisited),
    alpha=0.2
  )+
  geom_sf(
    aes(fill = factor(animal_ring_id)),
    data=spotted_nutcrackers_clustered_mcp
  ) +
  scale_animal_id() +
  scale_lv95_axes()


p1 + 
  p2 + 
  plot_layout(widths = c(5,5)) + 
  plot_annotation(tag_levels = "A")

Figure 3: Identified frequently revisited sites classfied by individuals as (A) points and (B) 95% minimum convex polygones

Figure 4 shows the same frequently revisited GPS locations, now grouped by animal sex. Both female and male individuals exhibit consistent use of distinct areas, particularly in the southeast and northeast of the study area.

p3 <- ggplot()+
  background_image(background)+
  geom_sf(
    aes(color=animal_sex), 
    data= spotted_nutcrackers_points|> 
      filter(!freq_revisited),
    alpha=0.2
  )+
  geom_sf(
    aes(color=animal_sex), 
    data=spotted_nutcrackers_points |> 
      filter(freq_revisited)
  )+
  scale_animal_sex()+
  scale_lv95_axes()+
  theme(legend.position = "none")

p4 <- ggplot() +
  background_image(background)+
  geom_sf(
    aes(color=animal_sex), 
    data=spotted_nutcrackers_points |> 
      filter(!freq_revisited),
    alpha=0.2
  )+
  geom_sf(
    aes(fill = factor(animal_sex)),
    data=spotted_nutcrackers_clustered_mcp
  ) +
  scale_animal_sex()+
  scale_lv95_axes()

p3 + 
  p4 + 
  plot_layout(widths = c(5,5)) + 
  plot_annotation(tag_levels = "A")

Figure 4: Identified frequently revisited sites classfied by animal sex as (A) points and (B) 95% minimum convex polygones

The analysis of variance indicates significant differences regarding stepwidth among individuals (p = 9.21e-11) as well as between sexes (p = 3.9e-05; Figure 5), with female animals showing a slightly greater stepwidth than male.

p5 <- spotted_nutcrackers_stepwidth |> 
  ggplot(aes(x = stepwidth, 
             y = animal_ring_id |> fct_rev(), 
             fill = animal_sex)) +
  geom_density_ridges(scale=2) + 
  geom_text(data = spotted_nutcrackers_aov_label, aes(x = 1, y = animal_ring_id, label = group),
            inherit.aes = FALSE, size = 4)+
  labs(
    y="density per individual"
  )+
  scale_x_log10(
    limits=c(1,20000),
    breaks=c(1,10,100,1000,10000)
  )+
  scale_fill_brewer(
    palette = "Spectral",
    name = "animal sex",
    labels = c("female","male")
  )+
  theme(
    legend.position = "none"
  )+ 
  expand_limits(y = length(unique(spotted_nutcrackers_stepwidth$animal_ring_id)) + 2)

p6 <-spotted_nutcrackers_stepwidth |> 
  ggplot(aes(animal_sex,stepwidth,fill=animal_sex))+
  geom_boxplot()+
  stat_compare_means(method = "anova", label.x.npc = 0.05)+
  scale_animal_sex()+
  scale_y_log10()+
  labs(
    x="animal sex"
  )

p5 + 
  p6 +
  plot_layout(widths = c(7.5,2.5)) + 
  plot_annotation(tag_levels = "A")

Figure 5: Differences in stepwidth between (A) individuals indicated by the density of their stepwidth and group according to the ANOVA and (B) sex

Discussion (simon)

References

Bracis, Chloe, Keith L. Bildstein, and Thomas Mueller. 2018. “Revisitation Analysis Uncovers Spatio-Temporal Patterns in Animal Movement Data.” Ecography. https://doi.org/10.1111/ecog.03618.

Calenge, Clement. 2024. adehabitatHR: Home Range Estimation. https://CRAN.R-project.org/package=adehabitatHR.

Graf, Valentin, Thomas Müller, Martin U. Grüebler, Urs G. Kormann, Jörg Albrecht, Anne G. Hertel, Marjorie C. Sorensen, Matthias Tschumi, and Eike Lena Neuschulz. 2024. “Individual Behaviour Shapes Patterns of Bird-Mediated Seed Dispersal.” Functional Ecology 38 (5): 1032–43. https://doi.org/10.1111/1365-2435.14556.

Maechler, Martin, Peter Rousseeuw, Anja Struyf, Mia Hubert, and Kurt Hornik. 2024. Cluster: Cluster Analysis Basics and Extensions. https://CRAN.R-project.org/package=cluster.

Sorensen, Marjorie C., Thomas Mueller, Isabel Donoso, Valentin Graf, Dominik Merges, Marco Vanoni, Wolfgang Fiedler, and Eike Lena Neuschulz. 2022. “Scatter-Hoarding Birds Disperse Seeds to Sites Unfavorable for Plant Regeneration.” Movement Ecology 10 (1): 38. https://doi.org/10.1186/s40462-022-00338-1.

Appendix

Wordcount

wordcountaddin:::text_stats()